%% basis_loc
function Phi = basis_loc(V, x, y, z, basis_type, derivative_type)
Jacobi = (V(2:4,:)-V(1,:))';
J_inv = inv(Jacobi);
P = [x, y, z];
P_hat = (J_inv*(P-V(1,:))')';
x_hat = P_hat(:,1);
y_hat = P_hat(:,2);
z_hat = P_hat(:,3);
if strcmp(derivative_type, "x")
    Phi = basis_ref(x_hat, y_hat, z_hat, basis_type, "x");
else
    dxi = basis_ref(x_hat, y_hat, z_hat, basis_type, "dx");
    deta = basis_ref(x_hat, y_hat, z_hat, basis_type, "dy");
    dzeta = basis_ref(x_hat, y_hat, z_hat, basis_type, "dz");
    if strcmp(derivative_type, "dx")
        Phi = J_inv(1,1)*dxi + J_inv(2,1)*deta + J_inv(3,1)*dzeta;
    elseif strcmp(derivative_type, "dy")
        Phi = J_inv(1,2)*dxi + J_inv(2,2)*deta + J_inv(3,2)*dzeta;
    elseif strcmp(derivative_type, "dz")
        Phi = J_inv(1,3)*dxi + J_inv(2,3)*deta + J_inv(3,3)*dzeta;
    else
        error("Invalid derivative type.");
    end
end
end

%% basis_ref
function N = basis_ref(x, y, z, basis_type, derivative_type)
if strcmp(basis_type, "P1")
    N = zeros(size(x,1), 4);
elseif strcmp(basis_type, "P2")
    N = zeros(size(x,1), 10);
else
    error("Invalid basis type.");
end
if strcmp(basis_type, "P1")
    if strcmp(derivative_type, "x")
        N(:,1) = 1-x-y-z;
        N(:,2) = x;
        N(:,3) = y;
        N(:,4) = z;
    elseif strcmp(derivative_type, "dx")
        N(:,1) = -1;
        N(:,2) = 1;
        N(:,3) = 0;
        N(:,4) = 0;
    elseif strcmp(derivative_type, "dy")
        N(:,1) = -1;
        N(:,2) = 0;
        N(:,3) = 1;
        N(:,4) = 0;
    elseif strcmp(derivative_type, "dz")
        N(:,1) = -1;
        N(:,2) = 0;
        N(:,3) = 0;
        N(:,4) = 1;
    end
elseif strcmp(basis_type, "P2")
    lambda = 1-x-y-z;
    if strcmp(derivative_type, "x")
        N(:,1) = lambda.*(2*lambda-1);
        N(:,2) = x.*(2*x-1);
        N(:,3) = y.*(2*y-1);
        N(:,4) = z.*(2*z-1);
        N(:,5) = 4*y.*z;
        N(:,6) = 4*x.*y;
        N(:,7) = 4*x.*z;
        N(:,8) = 4*y.*lambda;
        N(:,9) = 4*z.*lambda;
        N(:,10) = 4*x.*lambda;
    elseif strcmp(derivative_type, "dx")
        N(:,1) = 1-4*lambda;
        N(:,2) = 4*x-1;
        N(:,3) = 0;
        N(:,4) = 0;
        N(:,5) = 0;
        N(:,6) = 4*y;
        N(:,7) = 4*z;
        N(:,8) = -4*y;
        N(:,9) = -4*z;
        N(:,10) = 4*(lambda-x);
    elseif strcmp(derivative_type, "dy")
        N(:,1) = 1-4*lambda;
        N(:,2) = 0;
        N(:,3) = 4*y-1;
        N(:,4) = 0;
        N(:,5) = 4*z;
        N(:,6) = 4*x;
        N(:,7) = 0;
        N(:,8) = 4*(lambda-y);
        N(:,9) = -4*z;
        N(:,10) = -4*x;
    elseif strcmp(derivative_type, "dz")
        N(:,1) = 1-4*lambda;
        N(:,2) = 0;
        N(:,3) = 0;
        N(:,4) = 4*z-1;
        N(:,5) = 4*y;
        N(:,6) = 0;
        N(:,7) = 4*x;
        N(:,8) = -4*y;
        N(:,9) = 4*(lambda-z);
        N(:,10) = -4*x;
    end
end
end